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Abstract. Non-linear effects in the dynamical evolution of a shearing sheet made of stars are studied. First the 
implications of hitherto neglected non-linearities of the Boltzmann equation for the dynamical evolution of the 
shearing sheet are investigated. Using a formalism developed previously on the basis of the linearized Boltzmann 
equation it is demonstrated that the inclusion of the non-linear term leads to a feedback cycle for swing amplified 
density waves in the unbounded shearing sheet. Such a feedback is unique to star disks and is not known for 
gas disks. In order to present concrete examples of the non-linear feedback cycle a SCF code was developed and 
numerical simulations of the dynamical evolution of the shearing sheet are performed. The numerical results seem 
to confirm the theoretical predictions. The evolution of the shearing sheet resembles closely and might actually 
explain the recurrent spiral instabilities found in large-scale numerical simulations of the dynamical evolution of 
galactic disks. 
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1. Introduction 

The shearing sheet (Goldreich & Lynden-Bell 1965, Julian 
. & Toomre 1966) model has been developed as a tool to 
study the dynamics of galactic disks and is particularly 
well suited to theoretically describe the dynamical mech- 
anisms responsible for the formation of spiral arms. For 
simplicity, the model describes only the dynamics of a 
patch of a galactic disk. It is assumed to be infinitesimally 
thin and its radial size is assumed to be much smaller than 
the disk. Polar coordinates can be therefore converted to 
pseudo-Cartesian coordinates and the velocity field of the 
differential rotation of the disk can be approximated by a 
linear shear flow. These simplifications allow an analytical 
treatment of the problem, which helps to clarify the un- 
derlying physical processes operating in the disk. In the 
previous papers of this series (Fuchs 2001a, b, 2004, 2005) 
we have studied various aspects of the dynamics of spiral 
density waves in a shearing sheet of stars based on the lin- 
earized Boltzmann equation. In the present paper we study 
non-linear effects in the dynamical evolution of the shear- 
ing sheet. In particular, we investigate the role of the hith- 
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erto neglected non-linear term of the Boltzmann equa- 
tion. The principal result is that the non-linearity leads 
to a feedback cycle for swing amplification of spiral den- 
sity waves. Recurrent spiral instabilities, which lead to an 
ever changing appearance of the disk, have been reported 
by Sellwood & Carlberg (1984) and Sellwood (1989) from 
their numerical simulations of the dynamical evolution 
of self gravitating, differentially rotating disks. Toomre 
(1990) and Toomre & Kalnajs (1991) have argued that 
these recurrent spiral density waves are essentially due to 
swing amplified (Toomre 1981) random fluctuations of the 
surface density of the disks. Sellwood (1989), on the other 
hand, points out that contrary to expectation the spiral 
instabilities stay at the same level of activity in the exper- 
iments, if the particle number is increased and thus the 
particle noise is reduced. Moreover, Sellwood & Carlberg 
(1984) have shown that, once the cycles of recurrent tran- 
sient spiral density waves have established themselves in 
their simulations, the leading density wavelets which are 
then swing amplified are extremely unlikely due to chance 
alignments of the stars, because the amplitudes of these 
leading wavelets are too large. They consider the recur- 
rent spirals not to be unconnected random events aris- 
ing from swing amplified random noise. Sellwood (2000) 
puts forward the notion of a feedback mechanism for spi- 
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ral instabilities, which is related to fine structure of the 
distribution function of stars in phase space carved in by 
transient spiral density waves, which can then in turn in- 
cite density waves again. As a concrete example Sellwood 
& Lin (1989) and Sellwood & Kahn (1991) have devel- 
oped the concept of recurring groove instabilities. These 
rely on fairly sharp boundaries of the dynamically active 
parts of the self gravitating disks, because they are ini- 
tiated by edge modes (Toomre 1981). The groove insta- 
bilities develop then as rigidly rotating spiral modes of 
the disk 1 and seed fresh grooves through particle wave res- 
onance at their Lindblad resonances. However, Sellwood 
& Carlberg (1984) have shown in detail that the spiral 
structures seen in their numerical simulations are transient 
swing amplified shearing density waves exactly of the kind 
found in the infinite shearing sheet model. They demon- 
strate not only that the spiral arms shear from leading to 
trailing orientation, but also that the density waves grow 
preferentially with azimuthal wave numbers predicted by 
the shearing sheet model, which has become widely known 
as the Toomre (1981) "A — 2" prescription. Exactly the 
same was found by Fuchs & von Linden (1998) in their 
numerical simulations of the dynamical evolution of a self 
gravitating disk made of stars and interstellar gas. Toomre 
(1990) has argued convincingly that the spiral arms of Sc 
galaxies are such transient shearing density waves. One of 
the deeper reasons for this success of the infinite shear- 
ing sheet model in describing spiral density waves realis- 
tically is the rapid convergence of the Poisson integral in 
self gravitating disks (Julian & Toomre 1966). Consider, 
for example, the potential of a sinusoidal density pertur- 
bation 



/oo poo 
dx' / dy' 
-oo J —oo 



S 10 sin (kx') 



y/(x- x') 2 + (y-y' 



,(1) 



where G denotes the constant of gravitation. At y = 



$ = -4GE 10 sin(fc 2 ;) lim = 



27rGEio sin (kx) 



k 



(2) 



The sine integral in equation (2) converges so rapidly that 
at kxi, = 77 it reaches already 87% of its asymptotic value. 
Thus the "effective range" of gravity is only about a quar- 
ter of a wave length. The shearing sheet models effectively 
patches of galactic disks of such size. Thus the shearing 
sheet model seems to be a well suited tool to interpret 
the recurring transient spiral instabilities seen in the nu- 
merical simulations of self gravitating disks. Of course, the 
shearing sheet model with its rather drastic simplifications 
such as the neglected curvature of the mean guiding cen- 
ters of the stellar orbits and the assumed homogeneous 
mass distribution or the periodic boundary conditions in 



1 It can be shown (Toomre & Kalnajs 1980, unpublished re- 
sults, Sellwood & Kahn 1991, Fuchs 2005) that even the shear- 
ing sheet develops rigidly propagating, exponentially growing 
spiral modes, if a sharp left hand boundary is introduced. 



numerical simulations of the dynamical evolution has its 
limitations. Sellwood & Carlberg (1984) noted in partic- 
ular that the maximum growth factors of the shearing 
density waves in their simulations are somewhat larger 
then the predictions of the local model. Similarly Toomre 
(1981) finds also greater amplification of density waves in 
the full disk model of Zang. 

It is the aim of the present paper to demonstrate on 
the basis of the infinite shearing sheet model that a feed- 
back cycle of swing amplified density waves is indeed the- 
oretically expected. In the next section a formal deriva- 
tion is presented, while in the third section we describe 
numerical SCF simulations of the dynamical evolution of 
the shearing sheet and illustrate the feedback cycle with 
concrete examples. The SCF simulations are directly com- 
pared with N-body simulations of the dynamical evolution 
of the shearing sheet by Toomre (1990), Toomre & Kalnajs 
(1991) and Chiueh & Tseng (2000). 

Non-linear interaction of density waves has been pre- 
viously studied by Tagger et al. (1987) and Sygnet et 
al. (1988). These authors considered, however, the inter- 
action of rigidly rotating, slowly growing spiral modes, 
which do not develop in the unbounded shearing sheet. A 
preliminary report on the work presented here was given 
in Fuchs (1991). 



2. Theoretical considerations 

2.1. Boltzmann equation 

The time evolution of the distribution function of stars 
in phase space, /, is determined by the collisionlcss 
Boltzmann equation 



(3) 



where the square bracket indicates the usual Poisson 
bracket and H is the Hamiltonian of the stellar orbits 
(cf. paper I). Spiral arms are usually thought to be only 
minor perturbations of galactic disks. Thus we choose as 
in the previous studies a perturbation 'Ansatz' of the form 



f = fo + Sf, H = H Q + 5$. 



(4) 



where /o and Hq refer to the unperturbed shearing sheet, 
respectively. The physical idea is to subject the shear- 
ing sheet to a small spiral-like potential perturbation 5$ 
and to determine the disk response Sf by solving the 
Boltzmann equation. Inserting the 'Ansatz' (4) into the 
Boltzmann equation (3) leads to 



dt 



+ lf ,5<S>} + lSf,H }=-[5f,6$}, 



(5) 



where the rhs represents the non-linear term of the 
Boltzmann equation. The background distribution func- 
tion f is chosen in the following, again in the form of a 
Schwarzschild distribution function. 



B. Fuchs, C. Dettbarn, and T. Tsuchiya: Density waves in the shearing sheet 



3 



The non-linear Poisson bracket has the explicit form The linearized Boltzmann equation 



[5f,8*] = 

dw\ dJ\ dw 2 dJ 2 dJi dw\ dJ 2 dw 2 

if we use again the variables introduced in paper (I). J\ 
denotes the radial action integral of a stellar orbit and J 2 
is an integral of motion related to angular momentum. w\ 
and w 2 are the canonical conjugate variables, respectively. 
In expression (6) the perturbation of the distribution func- 
tion can be Fourier analyzed as 

Sf = J2 I dl 2 fi(J)e l[hwi+l2W2] . (7) 
h J 

Since w\ is a genuine angle variable, l\ denotes integer 
numbers, whereas l 2 is a real variable. The perturbation 
of the gravitational potential is Fourier analyzed as 



(5$ = J dk x J dk y ^e l[Kx+k y y] 



(8) 



where x and y denote the pseudo-Cartesian spatial coordi- 
nates of the shearing sheet with y pointing in the direction 
of rotation. Using the orbit equations derived in paper (I) 
the Fourier expansion in equation (8) takes the form 

J 2 — fio 



exp i [k x x + k y y] = cxp i (je x - 



/2Ji . 

+k x \j sin Wi + k y w 2 — /c. 



-2B 
2B 



■ COS W\ 



(9) 



with Oort's second constant B and the epicyclic frequency 
k. We refer to paper (I) for the meaning of f2o and ro, 
which define the center of the shearing sheet. Inserting 
the two Fourier components from equations (7) and (8) 
into the non-linear Poisson bracket (6) leads to 

[*/,<**] = 
., dfi .dfi \ 

- lk *dT 2 - l dTj cos{wi - w) r k 



xexpi^/iwi + l 2 w 2 + k 



J 2 - Oprp 
-2B 



+k y w 2 + £sin (wi - wj^J , (10) 

where the auxiliary variables £ and W are defined as in pa- 
per (I)as£ = \^\jkl + ^ky andW=arctan(^^), 
respectively. 



2.2. Solution of the Boltzmann equation 

The Boltzmann equation (5) can be viewed as a linear 
partial differential equation for Sf with the two inhomo- 
geneous terms [/o, <5$] and [Sf, <5$], respectively. Formally 
it can be solved for each inhomogcncity separately and 
both solutions can be then combined as the final solution 



Sf = Sfi + Sf 2 . 



(11) 



^ + [Sf 1 ,Ho} = -[fo,S^} 



(12) 



has already been solved in paper (I) and we can use di- 
rectly the result obtained there. This leaves the second 
equation 

dSf 2 



dt 



+ [Sf2,H 



-[5f,S$] 



(13) 



to be treated here, which takes the explicit form 



06 f 2 , dSf 2 



dt 



dwi 



5<* 



Q r ) 



dSf 2 
dw 2 



h jj- sin ( w i - w ).h + 



h 



( 

xexpi^hw-i + l 2 w 2 + fc x — — — — 
+k y w 2 + £ sin (wi - w)j 



(14) 



with the first Oort constant A = Qq + B. Equation (14) 
is of the same form as equation (20) of paper (I) and can 
be treated in the same way. Since the coefficients of the 
differential equation (14) do not depend on time t or the 
variable w 2 , equation (14) can be Fourier transformed with 
respect to t and w 2 leading to 



. , df 2oj 

iuUu + k— h i(ky 

dwi 



(-h - £Vo) /2,w = 



f £ fcx 



«TT7- < J , kU£cos (wi 
OJl 



/l*k|u> 

')) 



w) 



xexpi^/iwi + l 2 w 2 + fc x 



J 2 - Q ro 



-2B 



+k y w 2 + £ sin (wi - w)j . 



(15) 



Notice that the products of f\ and <f> k and -^j^ and <I> k , 
respectively, have been Fourier transformed with respect 
to time simultaneously, which has to be taken into ac- 
count when transforming back from frequency to the time 
domain. Equation (15) is of exactly the same form as equa- 
tion (23) of paper (I) and can be solved in the same way 
by determining first the solution of the homogeneous part 
of the equation, then constructing a particular integral 
of the inhomogeneous equation by 'variation of the con- 
stant', and finally adjusting the constant of integration so 
that the solution is uniquely defined in velocity space. We 
report here only the final result 

.l2.oj = ~^T- • expi[wt + (ky + l 2 )w 2 + k x — — 

2z sin 7177 — 25 

x (^fek/iSkL+i^lj-Sklu;) 

f° 

x / dw[ exp i[(l\ + r/)w[ + £ sin (w' 1 + w\ — w)] 

J-2-K 
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+h-z-rfl$k\w dw' 1 sin{w[+w 1 -w 

2J 1 J-27T 

x cxp +r])w' 1 + £ sin (w[ + Wi —w)] 
+i£^-j-$k| w / <M cos (w[ +wi-w) 

<J J l J-2-K 

x cxp + rj) w'x + £ sin (w^ + tui — to)] 



(16) 



where r\ = i + (fc y + Z2) 5 ( Ji — Slo T o) ) has been in- 
troduced as an abbreviation to keep the formulae a little 
more compact. w[ is an auxiliary variable. 

2.3. Self consistent perturbations 

Even in the linear case the spatial pattern of the disk 
response is different to that of the imprinted potential 
perturbation (cf. equations 10 of paper I). In particular, 
through the dependence on 77 and thus on the J 2 vari- 
able there is a further dependence of the disk response 
on the x coordinate. This means that the disk response 
to a Fourier component of the potential is not a Fourier 
component of the general disk response. Kalnajs (1971) 
has shown a way to overcome this difficulty by integrating 
the disk response over all wave numbers and then Fourier 
transforming this with respect to the spatial x and y coor- 
dinates again. The result is integrated over velocity space 
in order to obtain the surface density of the disk response 
and inserted into the Poisson equation whose lhs is also 
Fourier transformed, 



(27T) 



— 2 / dz dx 

7T ) J — 00 J — 00 J — 00 



dye 



-i[k' x x+k' y] 



Q2 Q2 Q2 r+00 , + oc 



x exp 



i[k x x + k y y] - ^Jkl + k^\z 

= -2^fc2 + fc 2 $k | w e i[k x x+k y y] 
(2tt) J-oo J -oo 



(27T) 

+ OC Z' + OO 

dk x I dk 

-00 



+ oc r + oo 

^ / du I dv5f . 

00 J —00 



(17) 



In this way self consistent perturbations of a self gravi- 
tating disk can be calculated (cf. equations 35 and 36 of 
paper I). Accordingly, on the rhs of equation (17) a multi- 
ple integral over the disk response to the non-linear term 
of the form 



AttG f + °° f + °° f 27T f^°° 
2 / dJi / dJ 2 I dw\ I dw2 

(2tt) Jo J-00 Jo J-00 

00 Z' + OC /* + 0O r- 

, i = _ 00 J-oo J-00 J- 



+ 00 



dk v 



exp 



-2B 



+ k' y W2 + £' sin (wi — w') 



x - . . expt[ujt + (ky + t 2 )w 2 + k x - 



2i sin7r?7 

1 dfx 



-2B 



x / dw[ cxp + f])w[ + £ sin (w[ + w\ — w)\ 
J -lit 



dwi sin (w[ + wi — w) 



x exp + r])w' 1 + £ sin (w[ + w\ — w)] 

— $kU / dw[ COS (w[ + W! -w) 
V J 1 J-2-K 

x cxpz[(/i + r\)w' x + £ sin (w^ + w± — w)] 



(18) 



has to be considered. Again we have changed from spatial 
(x, y) and (u, v) velocity coordinates to (Ji, J 2 ) and (wi, 
W2) coordinates. £' and and w' arc defined in an analogous 
way to £ and w, respectively, for primed wave numbers. 
The integrations with respect to J\ and J2 cannot be car- 
ried out directly, because /i(J) is not known cxplicitcly. 
Integrations over the angle variables w\ and w[ lead to 
intermediate results which arc not enlightening. But the 
integration over the W2 variable gives immediately 
f+00 

dw 2 e^ k y +h - k > 2 =5(k y + l 2 - fc y ) . (19) 



1 

2^ 



The integration over the wave number k y contracts the 
integrand of (18) to fc y = k' y — l 2 , which is the custom- 
ary sum rule of wave numbers well known from plasma 
physics. There is an analogous sum rule for the fc x wave 
numbers, which is, however, not directly obvious from ex- 
pression (18) because of the different Fourier transforms 
in equations (7) and (8). 



2.4. Volterra integral equation 

It was shown in paper (I) that the solution of the linearized 
Boltzmann equation (12) leads, upon inserting it into the 
Fourier transformed Poisson equation (17), to the Volterra 
integral equation (k' > 0) 



dkJZ (fc x , k' K ) $ k 



(20) 



where the kernel K. is given by equation (67) of paper (I). 
r l k , describes an inhomogeneity of equation (20) related 
either to an initial non-equilibrium state of the shearing 
sheet or an external perturbation. The formal solution 
of equation (13) as derived in the previous section and 
then inserted into the Fourier transformed Poisson equa- 
tion (17) has to be added to the rhs of equation (20) in 
order to obtain the solution of the non-linear Boltzmann 
equation. This means formally a further inhomogeneity of 
the Volterra equation (20), r^, which is given by (cf. equa- 
tion 36 of paper I) 



r k' — 



2G 



-2jk 



/• + OO p + OO /•2-7T 

/ dJi I dJ 2 I dwi 

Jo J-00 Jo 
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°° Z'+OO Z' + OO 

Y dl * dk * 

h = -oo J - x> " / -°° 

exp - i \k x — h 4 sin (wi - w'j 

x — — exp i [ut + k w 2 + fc x — — 

zism-Ki] y —2B 

X (^B ?2fcx/l$k|w+ * fcy |f $k| ") 

x / dw[ exp + r))w[ + £ sin (to^ + toi — to)] 

£ /-° _ 

+^i^-/i*kU / dtOi sin (to : + wi - to) 

x expi[(/i + rfjWi + £ sin (w[ + w\ — w)] 

+i£—— $ k | w / dw[ cos (to'j + ioi - w) 

OJl J-2-K 

x expi[(/i + ry)t«i + £sin (w[ + w\ — W)] . 



(21) 



The Volterra equation (20) is derived in the frequency 
domain. If transformed back to the time domain, we find 
for the linear part 

^,(t) = 2A\k' y \ f dt'K(k'^-2Ak' y {t-t')X) 

J — OO 

x$ ki -2Ak;(t-t'),k;(0 +rl(t) , (22) 

where the overhead tilde indicates that the w-dependent 
term of the kernel JC has been split off. In deriving equa- 
tion (22), use of the convolution theorem of the Fourier 
transform of products of functions has been made (cf. pa- 
per I) . In order to transform the non-linear term back into 
the time domain we consider first in expression (21) the 
terms that depend explicitely on ui, 

/oo+iv -i 
fajlvt+xv+Wi] _ _ = (23) 



sin (irn) 



-oo-\-iv 



oo+« — 



sin (tttj) 



with r] = i(cj + ^{J 2 — Qo r o)k' y )- According to Landau's 
rule the imaginary part of the frequency is chosen nega- 
tive, Ss(u>) = v < 0. By deforming the integration contour 
in the complex frequency plane we find that expression 
(23) is proportional to the unit step function 



u{nt + w[) 



1, Kt + w[ > 
0, Kt + w[ < ■ 



(24) 



Thus the inhomogeneity (21) of the Volterra equation can 
be cast formally into the form 



r-\-oc r-\-oo rZIZ 

r#(t) = 2G dJ 1 / dJ 2 / dwi 

JO J-oo JO 

Y dk * / dw 'i / dt> 

l 1 — _ OQ ^ — oo J —oo J —2ir J — OO 

[M(J, wi,w[,l, fcx, k', t - t')h(t')$ k ^ K - h (t') 



+AA 2 (J, w u w[, 1, fc x , k', t - t')^-\ t ,<t> KK - h (t') 

+Af 3 (J, Wl ,w[^k x y,t-t')^\t^K,K-h(t')] , (25) 

which has to be added to the rhs of equation (22). Again 
use of the convolution theorem of the Fourier transform 
of products of functions has been made. Since w[ is less 
than or equal to zero, the integration over time t' extends 
to the upper limit t. 



2.5. Discussion of the Volterra equation 

The Volterra equation (22) augmented by r£f reveals 
clearly the role of the non-linear term of the Boltzmann 
equation in the dynamical evolution of the shearing sheet. 

In Fig. 1 the linear evolution of a density wave, cal- 
culated neglecting the non-linear part of equation (22), is 
illustrated. The density wave is initialized as 



7v oc 6{k' x - (fci n + 2Ak' y t))S(k' y - k^e'^ , 



(26) 



where T is a time constant. As can be seen from Fig. 1 or 
directly from the interplay of the indices of the Fourier co- 
efficients of the gravitational potential in the linear equa- 
tion (22), k' K - 2Ak' y (t - t'), with time t' at which the 
Fourier coefficients are taken, the delta impulse travels 
at constant wave number k y n to positive wave numbers 
fc x at a speed of fc x ff = 2Ak y n . The wave crests of the 
density waves are oriented perpendicular to the wave vec- 
tors, which means that the wave crests swing around from 
leading to trailing orientation. While swinging around the 
wave is amplified and then dies out. In addition the ampli- 
tude shows oscillatory behaviour. Trailing density waves 
are not effectively amplified (cf. paper I). Since no outer 
boundaries have been introduced, no rigidly rotating spiral 
modes appear in the shearing sheet (cf. paper II). If the 
non-linear term is taken into account in equation (22), 
the evolution of the shearing sheet becomes more com- 
plicated. In the non-linear term (25) Fourier coefficients 
of the perturbations of the distribution function and the 
gravitational potential are combined. Since the non-linear 
term represents an inhomogeneity of the Volterra integral 
equation, each combination of the Fourier coefficients sam- 
pled over times t' earlier than t provides fresh input for the 
swing amplification mechanism. It is instructive to rewrite 
the inhomogeneity r£,' given by equation (25) formally as 



r#(t)= f dt S(t-t )r$(t 

J — OO 



(27) 



and to consider for each £o~epoch the corresponding 
Volterra equation 

<&u M {t)=2A\k' y \ f dt'jc(k' x -2Ak' y (t-t'),k' x ) 

J — OC 

(t') + 5(t-t )r^(t ). (28) 
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Fig. 1. Linear swing amplification of a density wave in 
the shearing sheet. The wave is initialized as a delta like 
impulse in wave number space with an initial wave vector 
(ft™, A;™) = ( — 2, 0.5) fc cr it, which corresponds to initially 
leading arms, and an initial amplitude of 0.04. The im- 
pulse travels at constant wave number k™ with an effec- 



tive radial wave number kt = k" 



2Ak] n t. The wave 



numbers are given in units of /c cr ; t defined below, and the 
time unit is S!^ 1 , the inverse of the mean angular velocity 
of the shearing sheet. Negative amplitudes at t > 9.5f2 ( ^~ 1 
are not shown. The parameters of the disk model are A/tto 
= 0.5 and Q = 1.4. 

Obviously the full Volterra equation is given by the super- 
position 

(t) = f dto*k',to(*) = 

J — oo 

2A\k' y \ I dt'jc(k' x -2Ak' y (t-t'),k' x ) 

x f dt ^ K _ 2Ak , {t _ t%Kito (t')+r^(t). (29) 

J — oo 

As can be seen from equation (28) each 5-impulse, sam- 
pling the history of the evolution of the sheet prior to io, 
is propagated linearly to time t independently of the fol- 
lowing evolution of the sheet. On the other hand, these 
linearly evolving 5-impulses influence later 5-impulses by 
non-linear interaction of the Fourier components, always 
cascading off fresh input for the swing amplification mech- 
anism. Only at to = t is there an instantaneous non- 
linear interaction of the Fourier components. Thus the 
non- linear interactions lead inter alia to a feedback cy- 
cle for swing amplification. Even the interaction of a pair 
of trailing Fourier components (k > 0), which are decay- 
ing, if coupled so that their wave vectors are subtracted 
from each other, leads to leading density waves as an in- 
put to the swing amplification mechanism, which is then 
amplified accordingly. Combined with some residual par- 
ticle noise this might well account for the recurrent spi- 
ral instabilities reported from the numerical simulations. 
As mentioned above, Sellwood & Carlberg (1984) have 
demonstrated by analyzing the evolution of individual spi- 
ral arms in their simulations that these are indeed swing 
amplified density waves of the same type as in the shearing 
sheet. The effect of the non-linear term of the Boltzmann 
equation thus has aspects that are completely different 
to the effects of non-linear terms in the hydrodynamical 



or Jeans equations. If the latter are used to describe the 
dynamics of density waves, the non linear advective and 
enthalpy terms imply instantaneous non linear interaction 
of the Fourier components and have been shown to lead to 
saturation and soliton-like self-modulation effects of the 
density waves (Ikeuchi & Nakamura 1976, Norman 1978, 
Meinel 1983). 

The feedback cycle described here is of course an effect 
depending quadratically on the magnitude of the pertur- 
bations of the sheet. Thus, if the amplitudes of the den- 
sity waves are at very low levels, the efficiency of the swing 
amplification mechanism might not be sufficient to sustain 
the cycle and it will die out eventually. 

3. Numerical self consistent field simulations 

In the previous section we have developed the theoret- 
ical concept of a feedback cycle for density waves in the 
shearing sheet mitigated by the non-linear coupling of the 
waves. In order to demonstrate concrete examples of this 
effect we have run numerical simulations of the dynamical 
evolution of the shearing sheet. 

3.1. Self consistent field method (SCF) 

The SCF method relies on a complete set of pairwise basis 
functions into which the density of the particles and the 
gravitational potential are simultaneously expanded. Each 
of the pairwise basis functions solve the Poisson equation. 
For the shearing sheet these are simply Fourier transforms, 

E(R) = y*d 2 fcS k e l ( k ' R ), (30) 

for the surface density of the sheet, and 

$(R,z = 0) = j d 2 k^e t( ^^ (31) 



with 

$k = — ;— £k 



(32) 



for the gravitational potential, where k = |k| (cf. equa- 
tions 32 to 34 of paper I). Changing to real quantities the 
transformation pair can be written as 

A(k) = J d 2 R E(R) cos (k, R) , A(-k) = A(k) , 

B(k) = J d 2 R E(R) sin (k, R) , B(-k) = -B(k) , 

1 f- OC p OO 

£(11) = —*/ dk x dfc y (A(k)cos(k,R) 

2tH J-oo Jo 

+B(k) sin (k,R)), 

/~( /'OO PC 

$(R, = 0) = / dk x 

n J-oo Ja 



dk v 



x - (A(k) cos (k, R) + B(k) sin (k, R)) . 

k 



(33) 
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Fig. 2. Initial set up of the 32768 particles in the simula- 
tion of the shearing sheet. The y-axis points into the direc- 
tion of the shear flow. The size is 2xq x 2y = 20 x 28A cr ; t . 



In the simulation the surface density is given by 



N 



j=i 



(34) 



where (xj,yj) denotes the position of particle j. AS is 
its contribution to the surface density, AS = Hq/N. 
Accordingly we find from equations (33) 



N 



A(k) = As cos ( k > R j) - 4S o 



sin (fe x i ) sin(fc y 7/ ) 



j=i 

N 



k x 



B{k) = As sin(k,Rj). 

j=i 



(35) 



In practice we cannot simulate an infinite sheet, but 
choose a finite size 2xo x 2y . This active sheet is sur- 
rounded by virtual sheets of the same size sliding along 
the central sheet according to the linear shear flow. If par- 
ticles are leaving the active sheet across one side, they 
are re-fed into the opposite side as if entering from one 
of the neighbouring virtual sheets with exactly the same 
particle distribution as in the active sheet (Wisdom & 
Tremaine 1988, Toomre 1990). Moreover, the SCF method 
is meant to simulate the fluctuations of the surface density 
of the sheet. Thus we have subtracted in equations (35) the 
Fourier coefficients of the unperturbed background density 



/xo /■yo 
dx / dy E e -*( k < R ) 
-x J -ya 
sin(k x x ) sin (k y y ) 



4S 



(36) 



Otherwise there would be an extra force pulling all parti- 
cles towards the center of the sheet. 

We have discussed in the introduction that the reason 
for the success of the infinite shearing sheet model describ- 
ing spiral density waves realistically is the rapid conver- 
gence of the Poisson integral in self gravitating disks and 
that the "effective range" of gravity of the density waves is 
only of the order A cr ;t/2. Thus, if the size of the simulated 
sheet is chosen large enough, it should behave like an infi- 
nite shearing sheet. Next the wave numbers are discretized 
with bin widths AA; X and Afc y , respectively. Discretizing 
equations (33) accordingly we obtain for the planar com- 
ponents of the gradient of the gravitational potential 

Gk n 



-Va>(Rj,z = 0) = ^Ak x ^Ak y 

kx,„ k y , n 

x [-A(k n ) sin (k n , Rj) + B(k n ) cos (k n , Rj)] , (37) 



n \l kx,n "f~ ^y,n 



with k n = (fc Xin , fc y ,n)- In eacn ti me step of the numerical 
algorithm the coefficients A(k n ) and B(k D ) and the forces 



(37) due to the fluctuations in the sheet are calculated 
anew. Then all particles are moved one time step ahead. 
The equations of motion are the epicyclic equations of 
motion augmented by the fluctuating forces (Toomre 1981, 
1990), 



x = 2Q y + 4£1 A - 



dx 



<9$ 

y = -2Q x - -jr(x,y) . 
dy 



(x, y), u = x, 
v = y + 2 Ax . 



(38) 



In a quiet disk the fluctuating forces would be zero and 
the particles would follow simply epicyclic orbits. In the 
simulations the equations of motion (38) are integrated 
numerically using the symplectic integrator described in 
the appendix. 

3.2. Trial runs 

The basic dynamics of the shearing sheet is characterized 
by Oort's constants A = B = A - Q , the 

epicyclic frequency of the particle orbits n = ^—4fl B, 
the epicyclic ratio of the velocity dispersions of the 
particles (Ju/ a v = 4fig/K 2 , the critical wave num- 
ber fccrit = k 2 /(2ttGT,q) and the critical wave length 
A cr it = 27r/fc cr it, and finally the Toomre stability param- 
eter Q = kct u /(3.36GS ). We introduce dimensionless 
model units in the following way. The length unit 
is A cr ;t = 1 and accordingly the wave number unit 
fc cr ; t = 2ir. We set GS = 1 which implies k — 2ir. The 
time unit is then the epicyclic period t K = 1. Moreover 
we assume a flat rotation curve, A = —B = ^Hq so 
that A = —B = and flo = \/2tt. The radial velocity 
dispersion is given by a u = ^fr<3 and the epicyclic ratio 
is a u /a v = V2. In all the simulations presented here we 
have used the following parameters: 
number of particles TV = 32768 
size #0=10, yo = V2x = 14.1 
time step 0.03125 

bin size of wave numbers Akx — Aky = 0.0625 
64 bins in fc x ranging from -4 to 4 
32 bins in fc y ranging from to 4 
Q=1.2 

The simulations were set up as illustrated in Fig. 2 ac- 
cording to Sellwood's (1983) quiet start recipe where the 
particles are distributed spatially homogeneously and ac- 
cording to a Schwarzschild distribution in velocity space. 
In Fig. 3 we show the simulation of the dynamical evolu- 
tion of the shearing sheet over 4 time units as snapshots 
of the 32768 particle positions. The sheet, which was ini- 
tially homogeneous, quickly develops density waves which 
can be seen as streaks in the frames in Fig. 3 in the time 
interval t ~ 1 to 2. However, the sheet then becomes fea- 
tureless again. This is due to rapid dynamical heating 
of the sheet which was already observed by Sellwood & 
Carlberg (1984) and in other numerical simulations of the 
dynamical evolution of galactic disks such as by Tomlcy 
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Fig. 3. Dynamical evolution of the shearing sheet. Snapshots of particle positions are shown for t = 0.5, 1, ... 3.5, and 
4, respectively. The t/-axis is oriented in the direction of the shear flow. The size of each frame is the same as in Fig. 2 



Fig. 5. Same as Fig. 3, but the shearing sheet is dynamically cooled by accretion of particles during the simulation. 



Of 




12 3 4 

time 



Fig. 4. The rise of Q' 2 = (u 2 )/Eq as function of time 
illustrating the disk heating. The time unit is the epicyclic 
period. 
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12 3 4 

time 

Fig. 6. Same as Fig. 4, but the shearing sheet is dynam- 
ically cooled by accretion of particles during the simula- 
tion. 



et al. (1991), Toomre (1990), Toomre & Kalnajs (1991), 
Chiuch & Tseng (2000). We illustrate the dynamical heat- 
ing in Fig. 4 where we show the evolution of the ratio 
Q' 2 = (u 2 ) /Sq, which mimics the Q parameter, as func- 
tion of time. During the evolution of the sheet its value 
has risen from 0.4 to 5.2, which reflects the rise of the ve- 
locity dispersion of the particles with time and that the 
sheet has become dynamically hot. The nearly linear rise 
of Q' 2 in the early phases of the simulation nicely fits the 
dynamical disk heating effect of density waves described 
theoretically by Fuchs (2001b) and Griv et al. (2002). In 
order to remedy this effect we have cooled the sheet dy- 
namically like Sellwood & Carlberg (1984) by adding dur- 
ing the simulation particles on orbits with their initial 
epicycle amplitudes. In Figs. 5 and 6 we show the results 
of a simulation which grew linearly from 16384 particles 
to 32768 particles as in the previous simulation. As can 
be seen from Figs. 5 and 6 the density wave activity and 
the pseudo-Q parameter stayed fairly constant through- 
out the simulation, even though the velocity dispersion of 
the particles was rising. Since the Q' parameter is larger 
than in the initial stages of the previous run, the ampli- 
tudes of the density waves shown in Fig. 5 arc smaller 
than in the early frames of Fig. 3. We illustrate the tips 
of the spiked amplitude distribution of the A(k) Fourier 
coefficients of the density waves in Fig. 7. The distribu- 



tion of the B(k) Fourier coefficients is very similar. As ex- 
pected from linear theory (Julian & Toomre 1965, Toomre 
1981, Fuchs 2001a) density waves with wave numbers in 
the region around k w (1, 0.5) & cr it implying pitch angles 
of about 27° have the largest amplitudes. These are the 
streaks seen in Fig. 5. They are, however, not station- 
ary patterns, but formed by the superposition of transient 
shearing density waves. A cut through the power spec- 
trum y A 2 (h) + B 2 (k) of the density wave amplitudes at 
a given k y wave number is also instructive. In Fig. 8 cuts 
at k y — 0.5 fc cr it are illustrated for the initial [t = 0) and 
a later stage (t = 2) of the simulation. Initially the power 
spectrum is flat in accordance with random distribution of 
the disk particles. Later on the power spectrum develops 
not only the high peak at positive fc x w 1 fc cr ;t wave num- 
bers which corresponds to the maximally swing amplified 
amplitudes of the shearing density waves, but is also signif- 
icantly higher at negative fc x wave numbers than the power 
spectrum of randomly distributed particles. Such leading 
density waves are the input to the swing amplification 
mechanism. We conclude in agreement with Sellwood & 
Carlberg (1984), who presented a similar argument, that 
the recurring density waves are indeed not amplified ran- 
dom noise, but are connected events. The enhanced power 
at leading wave numbers supports the argument that the 
coupling of trailing density waves can lead to leading input 
to the swing amplifier. 
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1 2 

K 

Fig. 7. Peaks of positive A(k) Fourier coefficients in the 
simulations presented in Fig. 5 at times t = 0.5, 1, ... 3.5, 
and 4, respectively. For clarity only contour levels of at 
least 40% of the maxima of the spikes are shown. The 
wave numbers are given in units of fc cr it- 

The SCF simulations shown in Fig. 5 look very similar 
to the direct N-body simulations of the dynamical evolu- 
tion of the shearing sheet by Toomre (1990) and Toomre 
& Kalnajs (1991). In these simulations an artificial fric- 
tion force acting on the particles in the sheet is included 
in order to overcome the disk heating problem. Toomre 
& Kalnajs (1991) explain the 'hotch-potch of swirling 
density waves' as the collective response of the shearing 
sheet to each of its individual members by polarization 
clouds around each particle. Instead of this quasi-linear 
concept we prefer to interpret our simulations in terms of 
the feedback cycle described in section (2). The SCF sim- 
ulations illustrated in Fig. 5 resemble closely the direct 
N-body simulations of Chiueh & Tseng (2000, their run 
a). Disk heating effects are apparently controlled there by 
the choice the softening length. 

3.3. Exciting a swing amplified density wave 

The numerically simulated shearing sheet is subjected to a 
"strong" external potential perturbation of the form con- 
sidered in section (2.5), 

(5$ k oc 6(k x - (fc x " + 2Ak y t))S(k y - jfej n )e~£ . (39) 

It was shown in paper (I) that according to linear theory 
the shearing sheet responds with a swing amplified shear- 
ing density wave. In Fig. 9 the result of a simulation with a 
time constant of T = 1.8 is traced in k-wave-vector space. 
As can be seen from Fig. 9 the disk response is indeed of 
exactly the same nature as predicted by linear theory. The 
impulse travels from negative fc x to positive wave numbers 
and the amplitude is maximally amplified around fc x s=a 2. 

3.4. Non-linear wave coupling 

The shearing sheet is now subjected to two independent 
impulsive external potential perturbations of the form 
given in equation (39) with initial wave vectors k*™' 1 and 
k m ' 2 , respectively. As shown in Figs. 10 and 11 both per- 




J i i i I i i i I i i i I i i i L 

-4 -2 2 4 

k x 

Fig. 8. Power spectrum y/A 2 (k x , 0.5) + 5 2 (fc x , 0.5) of the 
density wave amplitudes in the simulation shown in Fig. 5. 
The dashed and solid lines correspond to times t = and 
t = 2, respectively. The wave numbers are given in units 

Of fccrit. 

turbations evolve like the single perturbation in the pre- 
vious section. However, further density waves appear with 
circumferential wave numbers k y as predicted by the sum 
rule given in equation (19). We interpret these new density 
waves, which are delayed relative to the imposed perturba- 
tions, as examples of the onset of the non-linear feedback 
cycle predicted theoretically in section (2). In Fig. 10 an 
example is shown which started with the initial wave vec- 
tors k" 1 * 1 = (-2, 0.5) fc crit and k m ' 2 = (-3, 1) fc crit , respec- 
tively. At time t=0.88 there appears a new weak feature at 
about k = (0.9, 1.5) fc cr it- This then couples back with the 
first externally triggered wave leading to a strong wave at 
time t = 1.13 with k y = lfc cr iti which travels behind the 
second externally triggered wave. The example in Fig. 11 
started with initial wave vectors k 1 "' 1 = (—2, 0.5) fc cr it and 
k TO ' 2 = (—3, 1.5) fc cr ;t, respectively. At time £=0.5 a new 
weak feature appears at about k = (— 3, 1.5) fc cr it- After 
the second externally triggered wave was wrapped up very 
tightly, so that it travelled out of the frames shown in 
Fig. 11, there appears at time t = 1.13 a new strong wave 
at k = (1,1) fc cr it- Non-linear coupling of density waves 
with different wave numbers, which led to recurrent den- 
sity wave activity, has already been noted by Chiueh & 
Tseng (2000) in their simulations, although no detailed ex- 
planation was offered. They also find the decline of the am- 
plitudes of density waves at given fc y wave number which 
they ascribe to Landau damping. However, since density 
waves in the shearing sheet swing around with the shear 
flow, there is effectively no relative streaming of disk parti- 
cles relative to the waves (cf. Fig. 3 of Toomre 1981) , which 
would lead to Landau damping (Lynden-Bell & Kalnajs 
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Fig. 9. Response of the shearing sheet to an external impulsive potential perturbation with an initial wave vector 
(fc x ,fc y ) = (— 2, 0.5) fc cr it traced in wave number space. Frames are shown for time t = 0, 0.25,..., 1.5, and 1.75, 
respectively. 

Fig. 10. Response of the shearing sheet to two external impulsive potential perturbations with initial wave vectors 
(k x ,k y ) = (—2,0.5) fc cr it and (fc x ,fc y ) = (—3, 1) fc cr i t , respectively traced in wave number space. Frames are shown for 
time t = 0, 0.25, 0.5, 0.88, 1.0, 1.13, 1.25, and 1.38, respectively. 

Fig. 11. Same as Fig. 8, but with two external impulsive potential perturbations with initial wave vectors (fc x , fc y ) = 
(-2, 0.5) fc cr it and (fc x ,fc y ) = (— 3, 1.5) fc cr it, respectively. Frames are shown for time t = 0, 0.25, 0.5, 0.75, 1.0, 1.13, 
1.25, and 1.38, respectively. 



1972). Thus in our view hardly any Landau damping is 
to be expected. The density waves die out, because the 
density contrast is wiped out by the epicyclic motions of 
the stars, when the density wave crests are wrapped up 
tightly by the shear (Julian & Toomre 1966, Fuchs 2001a). 

4. Summary and conclusion 

The aim of this paper is to study non-linear effects in the 
dynamical evolution of an unbounded shearing sheet made 
of stars, which models a patch of a galactic disk. We have 
analyzed these effects both theoretically and by numerical 
simulations. 

First we have investigated theoretically the implication 
of the non-linear term of the Boltzmann equation for the 
dynamical evolution of the shearing sheet. The non-linear 
term is expressed in terms of Fourier expansions of the 
perturbations of the distribution function of the stars in 
phase space and the gravitational potential of the sheet, 
respectively, and the Boltzmann equation augmented by 
this term is formally solved. It is shown that this fur- 
ther term can be viewed formally as an inhomogeneity 
of the fundamental Voltcrra integral equation which de- 
scribes the dynamical evolution of the shearing sheet. The 
inhomogeneity is given by combinations of the Fourier co- 
efficients of the perturbations of the distribution function 
and the gravitational potential of the sheet, sampled over 
the past history of the evolution of the sheet. This leads 
to a feedback cycle for swing amplified density waves. It is 
quite different from the effect of non-linearities of the hy- 
drodynamical or Jeans equations, which have been used al- 
ternatively to describe the dynamics of the shearing sheet 
or galactic disks. 

In order to present concrete examples of the non-linear 
feedback cycle we have developed a self-consistent field 
code and have run numerical simulations of the dynamical 
evolution of the shearing sheet. The SCF method is based 
on a set of pairwise basis functions into which the sur- 
face density and the gravitational potential of the shearing 
sheet are expanded. Each of the pairwise basis functions 
solve the Poisson equation. In the case of the shearing 
sheet these are simply Fourier transforms. The equations 
of motion of the particles in the shearing sheet are the 
epicyclic equations of motion, which describe the motion 



of the particles in the unperturbed sheet, augmented by 
the forces due to the fluctuations of the sheet. In the simu- 
lations the shearing sheet turned out to be prone to rapid 
dynamical heating. Thus it was necessary to cool the sheet 
dynamically by adding mass continously to the sheet dur- 
ing the simulations. We have shown that the live shear- 
ing sheet responds to impulsive external potential per- 
turbations by developing shearing density waves. In par- 
ticular, by applying simultaneous multiple perturbations 
we have demonstrated the appearance of non-linearly in- 
duced shearing density waves, exactly as theoretically pre- 
dicted. 

The feedback cycle found here might well account for 
the hitherto unexplained recurrent shearing spiral insta- 
bilities seen in the large scale numerical simulations of the 
dynamical evolution of galactic disks. 
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Appendix A: A generalized leap-frog scheme 

We consider a set of ordinary differential equations 

ii = f(u,x), x = g(u,x), (A.l) 

where u and x may be also vectors. Integrating equations (A.l) 
by a small time step h from to to to + h leads to 

u(to + h) = M + uoh + itioft 2 + 0(h 3 ) 
= u + f(u, x)h + i/(u, x)h 2 + 0{h 3 ) 



+ y [g{uo, x ) ■ -7^) g(u , x ) + 0{K ) 



1 2 

= Mo + f(u,x)h + -h 
+0(h 3 ), 



'■*)'♦(»■£)' 



(A.2) 



a;(to + h) = xo + x h + -x h + 0(h ) 

— xo + g(u, x)h + ^g( u , x)h 2 + 0(h 3 ) 
1 2 

= xo + g(u,x)h + -h 
+0(h 3 ). 



( f --L) 9+ ( 9 -iL) 9 



(A.3) 



(<?(* 

= x + hg(u 1/2 ,x ) 

+ T {s{ui/ 2 , x ) ■ -^j g{u 1/2 ,x ) + 0(/i 3 ) , (A.6) 
and finally 

Ul/2 = U-l/2 + hf(u_ 1/2 ,X ) 

+ y (f(u-i/2,xo) ■ -^j f(u-i/z,x ) + 0{h 3 ) 
X! = xo + hg(u 1/2 ,x ) 

+ y (fl(«i/2, x ) ■ -^j g{u 1/2 ,x ) + 0(h 3 ) . (A.7) 

The equations of motion (36) of the particles in the shearing 
sheet can be written as a set of first order differential equations, 

m = 2D. v + / x 
v = 2Bm + f y 

x — u 

y = v-2Ax. (A. 8) 

From equations (A.7) we derive the leap-frog scheme 

Mi/2 = M-i/2 + h [2Q V- 1 /2 + fx(xo,yo)] 
+h 2 n [2Sm_ 1/2 +fy(x ,yo)] +0(h 3 ), 

Vl/2 = V-l/2 + h [2Sm_i/ 2 + fy(xo,Vo)] 

+h 2 (A-tlo) [2Q Q v_ 1/2 + f x (x ,yo)] + 0(h 3 ) , 

xi = xo + hui/2 + 0(h 3 ) 

yi = yo + h(v 1/2 - 2Ax ) - h 2 Au 1/2 + 0(h 3 ) . (A.9) 

The time step was set to h = 0.03125 in the simulations. Thus 
the orbits are integrated with an accuracy of h 2 = 0.001. 



Adopting the principle of the leap-frog algorithm (Hockney & 
Eastwood 1988) we introduce mid-points, 



X1/2 = x(t + h/2), Mi/2 = u(t + h/2) , 
and can write 



(A.4) 



/(mo, xi /2 ) = /(«o, xo) + - (g(uo,x ) ■ -^j f + 0(h 2 ) , 
g(u 1/2 ,xo) = 3 (m , xo) + ^ (j(uo,xo) ■ J^) g + 0(h 2 )(.A.5) 



It follows then 

M(t + h) = m + hf(uo,x 1/2 ) 

+ Y (/(«o,xo) • ^) f(uo,x ) + 0(h 3 ) 
= M + hf(u ,x 1/2 ) 

+ y (/("o, x 1/2 ) ■ /(m , X1/2) + 0{h 3 ) 
x(t + h) = xo + hg(u 1/2 ,x ) 
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